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Abstract 



We report on a Monte Carlo study of the lowest-Landau-level limit of the 
Lawrence-Doniach model for a layered superconductor. We have studied order 
parameter correlation functions for indications of the broken translational 
symmetry and the off-diagonal long range order present in the mean-field- 
theory vortex lattice. Our results are consistent with a single first order 
phase transition between a low temperature 3D vortex solid phase, with both 
broken translational symmetry and off-diagonal long range order, and a high 
temperature vortex liquid phase with no broken symmetries. We construct a 
universal phase diagram in terms of dimensionless parameters characterizing 
intra-layer and inter-layer couplings. The universal phase boundary extracted 
from our simulations and the associated latent heats and magnetization jumps 
are compared with experiment and with numerical results obtained for related 
models. 
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I. INTRODUCTION 



Because of the combination of high transition temperature, strong planar anisotropy, 
and short coherence length that occurs in high temperature superconductors, strong ther- 
mal fluctuations are present over a wide temperature interval in these materials. Thermal 
fluctuations are especially important in a magnetic field where they are responsible for the 
melting of the AbrikosovEl vortex lattice, often at temperatures that are well below the 
mean-field critical temperature, giving rise to a vortex liquid state^0 that is strongly dia- 
magnetic but also strongly resistive. In this paper0 we report on a Monte Carlo study of 
thermal fluctuations in layered superconductors with a magnetic field oriented perpendic- 
ular to the planes. Our study is based on a Lawrence-Doniach mo dell! and employs the 
lowest-Landau-level (LLL) approximationlli^El which is valid near the mean-field-transition 
temperature. The emphasis is on generic properties of this model rather than on quantitative 
estimates of phase diagrams, latent heats, and magnetization discontinuities for particular 
materials.0 We find that, for this model, a single first order phase transition occurs between 
a low-temperature state and a high temperature state with no broken symmetries. The low 
temperature state has off-diagonal long-range-order (ODLRO) both along and perpendicular 
to the field direction and broken translational symmetry. 

At the mean-field level, the phase transition between the normal state and the Abrikosov 
lattice state is highly unusual in two related respects. Firstly, the eigenvalue of the Cooper- 
pair density matrix which diverges at the transition point has a macroscopic degeneracy, 
in contrast to the isolated divergent eigenvalue with a zero-momentum-state eigenfunction 
found at zero magnetic field. The many divergent eigenvalues are associated with the many 
states in the Landau levels^ into which the transverse translational degrees of freedom of the 
Cooper pairs are quantized. At the mean-field level, this unusual property is responsible for 
the simultaneous development of superconducting order and broken translational symmetry 
in the Abrikosov vortex lattice state. Another consequence is that, at temperatures well 
above the mean-field transition temperature T C MF , fluctuations in a magnetic field in D- 
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dimensions are like those of a D — 2 dimensional syste at zero magnetic field. Brezin 
et a/.El first suggested, on the basis of expansions around the upper critical dimension for 
this problem (D = 6), that fluctuations would drive the Abrikosov transition first order. 
Indeed Hetzel et a/@ some time ago found evidence from Monte Carlo simulations that 
the transition is first order for D = 3. This conclusion has been substantiated by a large 
volume of subsequent!^ work from various points of view, although there is not yet universal 
agreement on the nature of the fluctuation altered phase diagram. In particular, some 
workers doubt the occurrence of any true phase transition in a magnetic field.il Others have 
suggested that the temperature T c at which superconductivity along the field first occurs 
is either largeiEflo or smallerBS than the temperature Tm where translational symmetry 
is broken and the vortex lattice forms. Here we find that, when the LLL approximation is 
valid, T c = Tm and the transition is first order. 

The LLL approximation is valid near the mean-field transition temperature. However, 
fluctuations in many high temperature superconductors are strong enough to drive the tran- 
sition field outside of its range of validity^ at least at the high temperatures where the 
transition field strengths can be reached in standard superconducting magnets. The LLL 
approximation is, by conservative estimates, accurate for fields exceeding half of the mean 
field H C 2. When fluctuations in higher Landau levels are still weak, their presence can be 
accounted for by renormalizing the parameters of the modelJiiE Partly for this reason, it has 
proved difficult to determine the range of validity of the LLL approximation by comparing 
with experimental data; for example the minimum field above which the LLL approximation 
is accurate in YBCO was estimated to be ~ 10 Tesla in one recent studyEl and ~ 2 Tesla in 
another.il A characteristic property of the LLL limit of the Lawrence-Doniach model is that 
the number of vortices passing through each plane is fixed by the magnetic field strength and 
does not fluctuate. The model must therefore fail qualitatively at very weak fields since, 
in the limit of zero magnetic field, thermally generated vortex loops are believed to play 
an essential role near the superfluid transition. We will return to this problematic issue in 
discussing the results of our simulations. 



We have previously reported on a simulation of LLL thermodynamics for two-dimensional 
electron systems.!! The present paper reports on simulations of the Josephson coupled lay- 
ers that constitute the Lawrence-Doniach model for superconductors with strong planar 
anisotropy. In Section II we summarize and discuss the LLL limit of the Lawrence-Doniach 
model and introduce the correlation functions whose temperature dependence we have stud- 
ied. In Section III we present our simulation results. In Section IV we discuss some im- 
plications of our simulations, commenting on their relationship to experiment and to other 
simulations. Finally we briefly summarize our study in Section V. 

II. LLL LIMIT OF THE LAWRENCE-DONIACH MODEL 

In the Lawrence-Doniach model Jll the local superconducting order parameter \I/ n (f) is de- 
fined on discretely labeled continuum layers. (Here and in what follows r is a two-dimensional 
(2D) coordinate, taken to be the x — y plane, and n is a layer index.) The Ginzburg-Landau 
free energy for a particular configuration of the order parameter is given by 

F GL = d J2 fd 2 ?mn] (1) 

n J 

where 

im = £HT)|* n (f)| 2 + § |* n (r)| 4 + -^-\(-ihV - 2eA)* n (r)\ 2 + 
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■-— Hl*n + i(r)exp(— A z dz)-* n (r)\% (2) 
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do is the layer thickness, and d is the distance between layers. In the limit where typical 
order parameter functions vary smoothly as a function of layer index this reduces to a 
three-dimensional (3-D) Ginzburg-Landau model with anisotropyS parameter •y 2 = m c /m a b. 
The mean field transition temperature in the absence of a field (T^ F ) occurs where a(T) 
vanishes. Where definiteness is required we will follow the usual practice of assuming that 
a varies linearly with temperature and that other parameters of the model are temperature 
independent. (This assumption is unlikely to be valid across the entire fluctuation regime 



in typical high-temperature superconductors.) We limit our attention to magnetic fields 
directed perpendicular to the layers and choose a Landau gauge with vector potential A = 
(0, Hx, 0) where H is the magnetic field strength.0 The 2D gradient term in this energy 
functional is diagonalized by kinetic energy eigenfunctions for particles of charge 2e in a 
perpendicular magnetic field. In the underlying fermionic description these particles are the 
electronic Cooper pairs. Sufficiently close to the mean-field critical temperature, fluctuations 
are dominated by those in the Hilbert subspace of macroscopic dimension within which 
the Cooper pair kinetic energy is minimized. Within this manifold, the 2D kinetic energy 
operator can be replaced by a constant, heH/m a bC. This replacement, combined with the 
corresponding constraint on the order parameter, is the LLL approximation. 

Our simulations are performed on finite systems consisting of N z rectangular layers with 
sides of length L x and L y . The quasiperiodic boundary conditions consistent with our gauge 
choice are: 

^ n {x,y + L y ) = ^ n {x,y) 

^ n (x + L x , y) = exp(z2vr^(|-))^ n (x, y) 

Ly 

V n+Nz {x,y) = y n {x,y) (3) 

Here = (L x L y )/(2ir£ 2 ) is the number of flux quanta of the internal field passing through 
each layer which must be an integer if the finite-size boundary conditions are to be satisfied. 
(The magnetic length i 2 = hc/2eH.) In the LLL approximation, the order parameter in 
each layer is expanded in terms of a complete set of minimal-gradient-energy eigenfunctions 
that satisfy these boundary conditions. It turns out that this 2D basis set consists of 
displaced elliptic theta functions.il'il. Therefore the order parameter is expanded in the 
form 

*»(f) = ( l -^-) 1/2 J2CnA £ {Ly)-^{^)-^eM^yx 3 Ji 2 )eM-{x-x h ^ 2 /2i 2 )} 

j S = — 00 

(4) 

In Eq. © a H = a(T) + (heH/m ab c) = a(T)(l-H/H% F (T)), X jtS = j2ir£ 2 / L y + sL x , j runs 
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from 1 to iV^, and n runs from 1 to N z . The normalization is chosen so that the average value 
of \cj\ 2 is ~ 1 for Abrikosov's mean-field order parameter. As we see explicitly in Eq. (|]), 
the order parameter in the LLL approximation in this gauge is the product of exp(— x 2 /2£ 2 ) 
and an analytic function of z = x + iy. This constraint, whose analog in another gauge has 
been emphasized in the work of Tesanovic and collaborators^ and powerfully employed in 
studies of the quantum Hall effect, means that amplitude and phase fluctuations of the order 
parameter are necessarily linked in this regime. Since the line integral of the phase gradient 
of an analytic complex-valued function around any closed loop is 2ir times the number of 
enclosed zeroes, it follows from Eq. (H) that for each layer index n and each configuration of 
the order parameter, \l/ n (r) has precisely zeroes. There are no thermally generated flux 
loops in the LLL approximation. Moreover, for each n, ^ n {r) is specified up to a complex 
overall factor by the positions of its zeroes. 

We will characterize the states that appear in our simulations in terms of intra-layer and 
inter-layer correlation functions constructed in terms of the quantity A n / n (g) defined by the 
following equation: 

j d 2 r^> n ,(r)m n (r) exp{-iq ■ f} = A^A <w (g) exp{-^-}. (5) 

In terms of order parameter expansion coefficients for the finite system A n i <n (q) has the 
form: 

A n',n(g) = jfY, Cn'j 1 C n j a 6 ja -. jl - ny ex^[-iq x {X h + X h )/2] (6) 

3132 

where Sj = 1 if j is a multiple of and is zero otherwise, and Xj = Xj^- For the Fourier 
expansions of our finite systems with periodic boundary conditions, q = 2n(n x /L x , n y /L y ) 
where n x and n y vary from over any range of consecutive values. This quantity is 
conveniently sampled in our Landau gauge Monte Carlo simulations. A n / ( „(g) is proportional 
to the Fourier transform of an order parameter product that is diagonal in planar coordinates 
but off-diagonal in layer indices. We remark that because of the analyticity property of 
LLL wavefunctions this quantity nevertheless completely specifies the off-diagonal order 
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parameter product. In particular, inverting the Fourier transform and using analyticity it 
can be shown thatHMl 

= exp (-l^g) eM i{xy ~ £2 x ' y,) ) E «p(-^) 

x exp(^^ ^ exp(^^ 4. (7) 

We will principally be interested in whether or not our simulations show evidence for 
broken translational symmetry in the spatial distribution of the local superfluid density 
|^n(^)| 2 or evidence for ODLRO in the order parameter. We define the quantity 

A = ^£A„, n (g = 0). (8) 

Ao is proportional to the average local superfluid density over all layers. In mean-field theory 
A may be taken as the order parameter of the Abrikosov state. (With our normalization, 
A = 2//3 A wherel! /3a ~ 1.159595 in the Abrikosov state.) When fluctuations are included, 
we expect that A will vary from zero well above the mean-field transition temperature to 
2/ /3a at zero temperature, with singular behavior at an intermediate temperature if a phase 
transition occurs. A non-zero value for Ao does not indicate either broken translational 
symmetry or ODLRO. Broken translational symmetry can be identified by examining 

_ A <n /(-g)A n , n (g) 

Sn'AQ) = X2 • ( 9 ) 

^0 

S n t n measures intra-layer (n' = n ) and inter-layer (n' 7^ n ) superfluid density correlations. 
In the Abrikosov mean-field solution S n >,n{o) — 1 when q equals a reciprocal lattice vector 
(G) of the vortex lattice and is otherwise zero. For states with broken translation symmetry, 
the thermal average of S n ^ n {G) will remain finite for — > 00. For a vortex liquid state 
we expect that S n ^ n (q) will vanish as iV^ 1 at large A^, for all wavevectors q. S n ^ n (q) is 
tantamount to the static structure factor of the vortices, although this equivalence is not in 
general precise as far as we are aware. 

Off-diagonal long-range order along the field can be identified by examining 
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1 v A n> Qf=0) 

For the Abrikosov mean- field solution z — 1. For states with ODLRO along the field the 
thermal average of z will remain finite in the thermodynamic limit, otherwise it should go 
to zero as N~ . Application of the notion of ODLRO to directions perpendicular to the field 
is more subtle and will be discussed later. 

The Ginzburg-Landau free energy can be expressed in terms of A n / n (g*): 



1 Q 2 1 2 

sgn(a H )A n , n (0) + - £ |A n , n (q)| 2 exp(-^-) 
q 



n 

+2r/(A„ in (0)-Re(A nin+1 (0)))] 
= N^N z g 2 (E 2d + V Ej). (11) 



Here g = any irPdo/ ^ksT and r\ = (h 2 /2m c d 2 )\aH\ 1 measure the intra-layer and inter- 
layer coupling respectively, E 2d = N' 1 En[sgn(a H ) A n , n (0) + Eq(| A n , n (q)| V 4 ) exp(-q 2 £ 2 /2)] 
is the intra-layer energy and Ej = 2N~ 1 r]J2 n (^n,n(^) — Re(A njn+1 (0))) is the inter-layer 
Josephson coupling energy. For an{T) < the square of the dimensionless intra-layer 
coupling constant (g 2 ) is the ratio of the mean-field-theory superconducting condensation 
energy per-vortex per-layer to the thermal energy ksT. For typical high temperature super- 
conductors this ratio remains small over a wide range of temperature below the mean field 
transition temperature. The dimensionless Josephson coupling constant r\ is the ratio of the 
square of the Gaussian approximation to the correlation length along the field direction, 

,2 = ^ = ,2 W(T)\ 

?H 2m c \a H (T)\ ?c K(T)j' 1 } 

to the square of the layer separation. Several researchersil'0 have speculated on possible ex- 
perimental consequences of the dimensional crossover from two to three dimensions, naively 
expected when 77 exceeds 1. We will comment on these issues later. 
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III. NUMERICAL MONTE CARLO SIMULATIONS 



We performed Monte Carlo simulations of the LLL limit of the Lawrence-Doniach model 
for finite size systems. The complex expansion coefficients C n j of the order parameter in 
Eq. dJ) were treated as classical variables in the simulations. The real and imaginary parts 
of C n j were independently incremented by random numbers from a distribution chosen so 
that half of all attempted moves were accepted using the standard Metropolis algorithm. 
Each Monte Carlo step consisted of update attempts for the real and imaginary parts of all 
N^N Z coefficients. We calculated distribution functions and averages for various quantities 
of interest, including the total Ginzburg-Landau energy, the separate planar and Josephson 
coupling contributions to the Ginzburg-Landau energy, interlayer and intra-layer local su- 
perfluid density correlation functions, and correlation functions that are off-diagonal in the 
layer index. The simulations were performed for a range of values of the two dimensionless 
parameters that characterize the system, g and rj, and for a range of system sizes. We have 
focused our attention primarily on simulations at small values of rj, where the discreteness of 
the superconducting layers plays an important role, because this is the situation of relevance 
to most high-temperature superconductors. 

Most of our Monte Carlo simulations were performed for finite systems containing N z = 8 
layers and = 16, = 24, or = 36, vortices per layer. The number of Monte Carlo 
steps used in a single simulation was typically M = 2 x 10 6 . Our finite size system shapes 
have been chosen to accommodate perfect triangular lattices by taking 

r = 7f x F (13) 

where N y and N x are integers and N$ = N x x N y . (We choose N x = N y = 4 for Nj, = 16; 
N x = 6, N y — 4 for = 24 and N x = N y = 6 for = 36). For all simulations the order 
parameter was initialized to the Abrikosov lattice value and the first 10 4 Monte Carlo steps 
were discarded. We have compared results from different stages of our Monte Carlo runs to 
make sure our systems were well equilibrated. 
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Some typical results for (S n / >n (q)) are shown in Fig. |]. All the results reported here 
are for temperatures below the mean-field transition temperature so that g < and larger 
\g\ corresponds to lower temperatures. In isolated 2D systems, the pancake vortices of 
the LLL model meltH at g ~ —6.6. We expect coupling between layers to increase the 
transition temperature. The results in Fig. [I| show that, for this system size, a qualitative 
change in local superfluid density correlations occurs between the two temperatures for 
which results are shown. At the higher temperature (g = — vTo), the peak in S n>n (q) near 
reciprocal inter- vortex distances demonstrates that the pancake vortices are in a moderately 
correlated liquid state. There is no evidence of correlations between vortices in different 
layers, perhaps not surprisingly in view of the small value for 77 = 0.01. On the other hand 
for g = — V30, S n / >n (q) is nearly independent of the layer indices, indicating that vortices 
in different layers are highly correlated. Moreover, S n > >n (q), is strongly peaked at reciprocal 
lattice vectors, reaching a value that is approximately 70% of the Abrikosov state value 
for the nearest neighbor shell. Apparently, once the vortices are strongly correlated within 
a layer, weak interlayer coupling is sufficient to strongly favor configurations in which the 
vortex coordinates are weakly dependent on layer index, i.e. configurations with nearly 
straight vortex lines. 

Our simulations indicate that the transition between these two patterns of correlation 
happens abruptly at what appears to be a first order phase transition, which for rj = 0.01 



occurs for g »s — V20.5. In Fig. |2| we show the probability distribution function for the 



Ginzburg-Landau free energy at r\ — 0.01 and g = — V20~5. The double-peaked distribution 
function with a intermediate minimum which deepens with system size is strong evidence!!!! 
that a first order phase transition occurs in the thermodynamic limit in this system. The 
logarithm of the ratio of the maxima in P(E) to the intervening minimum may be interpreted 
as a free energy barrier between the low temperature vortex lattice state and the high 
temperature vortex liquid state. The transition will be first order if this barrier diverges in 
the thermodynamic limit. By this measure, even for this relatively weak Josephson coupling, 
the phase transition is much more strongly first order than for uncoupled layer systems. In 
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that case the free energy barrier grows more slowly with system sizes and first appears in 
the simulations only when the number of vortices per layer is close to lOoHH. The energies 
in Fig. are given in units of the mean-field-theory condensation energy. We see that for 
this value of rj, the Ginzburg-Landau free-energy at the transition is, in magnitude, about 
93.5% of its mean-field value when in the vortex lattice state and about 90.5% of its mean 
field value when in the liquid state. It would be tempting to identify the separation between 
the two peak energies in Fig. |^ with the latent heat. However, as we discuss later, this 
identification would be incorrect because of the temperature dependence of the parameters 
of the Ginzburg-Landau free energy. 

Fig. ^contains a contour plot of the two-dimensional distribution function P(E 2 d, Ej) at 
the transition. This picture reinforces the picture of the transition as being controlled by a 
competition between planar and Josephson coupling energies, both of which favor order, and 
entropic contributions to the free energy which favor disorder. The planar energy decreases 
because typical vortex configurations more closely approximate their triangular lattice op- 
timum while Josephson coupling energies drop because of increased interlayer correlation. 
The distribution function plot demonstrates that both the planar Ginzburg-Landau energy 
and the Josephson coupling energy are lower in the ordered state. We have defined Ej so 
that it is zero if the order parameter is identical in every layer. For completely uncorrelated 
order parameters in different layers, r = Ej/A = 2. We find that F has the values 0.598 
and 1.208 in ordered and disordered states respectively at the r] = 0.01 phase transition; 
T drops from approximately 60% to 30% of its uncorrelated layer value on going form dis- 
ordered to ordered states. This large change should be compared with the more modest 
relative change in the planar condensation energy. We show later that the latent heat and 
magnetization discontinuity at the first order phase transition are related to the changes 
in planar condensation energy and Josephson coupling energy that can be read off these 
figures. 
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IV. DISCUSSION AND COMPARISON WITH EXPERIMENT 



A. Phase diagram 

We have completed Monte-Carlo simulations of our model at a series of values for the 
two dimensionless parameters (g and n) which characterize it. These two parameters can be 
related to superconductor material properties that are measurable, at least in principle: 



2 _ A T (T e0 ) T c0 T c0 H^ F (TV H 

A e// T [ T ' H 1 H^ F {T) 



and 



where A T (T) = $g/167r 2 fc B T = 2 x 10 8 A/T is the thermal fluctuation length, A e// = 
2A 2 b (0)/(io is the two-dimensional (2D) effective penetration depth, k = A/£, and 7 2 = 
m c /m a b. The numerical form of the equation for gn applies for H in units of Tesla, T in 
units of Kelvin, and do and d in units of A. Note that gn does not vary strongly on crossing 
through the fluctuation regime with varying temperature or field, and that this product can 
be quite small for strongly anisotropic extreme type-II superconductors. Where our model 
is valid, the phase diagrams and appropriately scaled physical properties (see below for 
examples) of all planar superconductors should become identical when field and temperature 
are expressed in terms of g and n. 

Fig. § shows the universal (n, g) phase diagram we have constructed from Monte Carlo 
simulation data like that discussed in the previous section. The construction of the phase 
boundary line, n m = f(g m ), is aided by the following Clapeyron^l identity: 

dn m 2 6 < E 2d > , , , 1fi v 

-j — = — h F ^Vml- (16) 

dg m g l 6 <Ej> 

Here 5 < E 2 d > and 5 < Ej > are the respective discontinuities across the phase boundary 
of the planar condensation energy and the Josephson coupling energy discussed in the previ- 
ous section. We have used points on the boundary and these derivatives to obtain the spline 
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fit to the phase boundary line that is indicated in Fig. |]. The data points used for this fit 
were (r)m,g m ) = (0.0,-6.6), (0.005,-5.2), (0.01,-4.5), (0.02,-3.9) and dr] m /dg m = 0.024 
at (r] m , g m ) = (0.02,-3.9). The dimensionless parameters of a particular layered super- 
conductor will move through this phase diagram from upper left to lower right on crossing 
through the fluctuation regime by decreasing temperature or field. The three additional 
curves in Fig. |4| are lines of constant grj appropriate for experiments performed at typical 
fields and temperatures on various materials. For the quasi two-dimensional BSCCO family 
of superconductors the typical value of rj at the phase boundary is quite small. (At H = 1 
Tesla, T = 50 Kelvin, r] = 0.0014 for BSCCO.) On the other hand, for bulk superconductors 
grj at the phase boundary can easily be larger than 1.8, the value of the product at the upper 
right of the region illustrated in Fig. |j. 

The two solid dots in Fig. [| at rj = 0.06 and r] = 0.10 were not included in determining 
the spline fit but lie quite close to the line extrapolated from weaker couplings. At these 
values of rj, our N z = 8 simulations found that the measured correlation functions changed 
abruptly with g, but did not uncover the double-peak distribution in the Ginzburg-Landau 
free energy that is characteristic of a first order phase transition. At 77 = 0.06 more layers 
are required before the simulation begins to show a double-peak distribution; in Fig. [S] we 
show results obtained at the phase boundary for N z = 32 and = 16. 

As 77 increases, the correlation length along the field direction increases, and finite N z 
effects should increase in importance. At larger values of 77, we believe that the phase 
boundary line would move toward the right if we increased the number of layers. One way to 
judge the importance of finite size effects is to recognize that for 77 — > 00 the order parameter 
will become layer independent for any finite N z and our model will reduce to a 2D model 
with g 2 — > g 2 N z . Since for an isolated 2D layer we have foundSS that g m (j) = 0) = —6.6, we 
expect to find g m = —Q.Q/\fW z whenever the correlation length along the field direction in 
the fluid phase has reached the system size N z d. For example, for simulations with N z = 8, 
like those on which Fig. |] is based, g m will never be less than rs 2.3. The phase boundary 
we find at 77 = 0.12 is perilously close to this limit and it is likely that for N z = 8 we are 
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quite close to the 2D limit at this interlayer coupling.El Even for YBCO it seems likely that 
more than N z = 8 layers are required for a simulation of the phase transition^. 

As mentioned previously, for materials with reasonably narrow fluctuation regimes, — grj 
is essentially constant on sweeping through the transition regime with either field or tem- 
perature. Values of —grj for a number of planar superconductors at typical fields and tem- 
peratures are indicated in Table |. With this in mind we remark that in high temperature 
superconductors, and in other highly anisotropic planar systems, the transition to the vor- 
tex lattice state will always occur in the regime where rj < 1 and the discrete nature of 
the superconducting layers is of importance. To see this we observe that in the limit of 
large rj our model reduces to the LLL model for continuum three dimensional systems. In 
that limit, properties of the model are, as in the isolated layer limit, controlled by a single 
dimensionless paramete r§ 

Ad = 9 2 VV (17) 

g\ D is the ratio of the mean-field-theory condensation energy, per vortex and per correlation 
length, to the thermal energy ksT. It follows from Eq. (|I7|) that for large rj, rj m oc g~ 4 
in Fig. |4j. The phase transition to the Abrikosov state must still be strongly influenced by 
fluctuations for rj ^> 1 since dimension D = 3 is still below the naive lower critical dimension 
for this problem. There is, as yet, no quantitative estimate of the value of g^o at which the 
putative phase transition occurs in this 3D limit, but it seems clear that its value must at 
least be larger than unity. If this is so, the value of rj at the phase boundary must satisfy 
the inequality rj m < (grj)^ 3 . Only materials with grj > 1 can be in the 3D continuum limit 
at the phase boundary. Our attention here is strictly focused on materials with grj < 1. 



B. Thermodynamics 

It is useful to discuss the thermodynamics of the LLL Lawrence-Doniach model using 
a microcanonical ensemble approach that we introduced earlier for isolated layersBS and 
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we briefly generalize below. This approach is based on the observation that the Ginzburg 
Landau free energy can be expressed in terms of a small number of physically meaningful, 
intensive thermodynamic variables. For the present case we introduce three variables: Aq, 
the Abrikosov ratio Pa 

fiA=N ** M) &j*wm* (18) 

and a Josephson coupling parameter 

r _ En/^r1vl/ rt+1 (f)-vl/ n (r)| 2 

The Abrikosov ratio is small for smooth order parameters whose zeroes are relatively evenly 
spaced in each plane; it's minimum value Pa = Pa is achieved when the vortices are placed 
on a triangular lattice!^ and it reaches the maximum value Pa = 2 for the vortex gas state, 
r has minimum and maximum values and 2 for the Abrikosov lattice and vortex gas states 
respectively. Fluctuations in these intensive thermodynamic variables will vanish in the 
thermodynamic limit. In Fig. |6|, we plot A , Pa, and T as a function g for a finite system 
with fixed —gt] = 0.1. The melting transition for this value of gt] occurs near g = —3.6 
where the Josephson coupling parameter decreases relatively abruptly. 

With these definitions the LLL Lawrence-Doniach Ginzburg-Landau free energy can be 
expressed in the form: 

F GL = k B TN+N g f (20) 

where 

/ = g 2 (sgn(g)A + ^ + ^AqT). (21) 

The free energy of the LLL Lawrence-Doniach model can be expressed as the sum of Fql 
and an entropic contribution proportional to the log of the volume in the model's phase 
space with given values of the three variables. Noting that both Pa and V have been defined 
so that they are invariant under scale changes of the order parameter, it follows^ that the 
free energy of the model is 
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F = k B TN^N z [f - ln(Ao) - s{(3 A , V)}. (22) 

where s(/3a, T) is an entropy function that completely specifies the thermodynamics of the 
model. The equilibrium values of the thermodynamic variables Ao, /3a, and T and the free 
energy are determined at specified values of the dimensionless model parameters, g and rj, 
by minimizing the right hand side of Eq. (|22|). For example, minimizing with respect to A , 
we find that for g < 0, i.e., below the mean-field transition temperature, 

2(1 -rjT) 2 

^ - — sr~ + w m 

This constraint on the three thermodynamic variables we have introduced will be used below 
to simplify expressions for the latent heat and the magnetization discontinuity at the model's 
first order phase transition. Minimizing with respect to T or (3a yields a relationship that 
involves the entropy function. It is easy to show that s(/3, V) reaches its maximum value at 
the high temperature equilibrium values (3a = 2 and T = 2; however the properties of this 
function have not yet been thoroughly investigated. Minimizing with the entropic terms in 
Eq. ( [22]) ignored produces mean field theory. 

The thermodynamic signature of a first order phase transition between the disordered 
vortex liquid phase and the ordered vortex solid phase is the presence of a latent heat and 
a discontinuity in the magnetization of the system. The latent heat is proportional to the 
decrease of entropy on entering the ordered state. In our LLL Lawrence-Doniach model the 
entropy can be expressed in the following form 

Since the free energy is continuous at the transition, the latent heat has contributions from 
the discontinuity in (Fgl) and the discontinuity in the last term on the right hand side 
of Eq. (0). The first contribution to the latent heat is what would have been read off the 
separation between the two-peaks in the distribution function of Fgl near the transition and 
represents a contribution from the order parameter phase space volume. The second term is 
due to the temperature dependence of Fgl and represents a contribution from fluctuations 
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on a microscopic level for a given order parameter function. This contribution to the entropy 
is proportional to the magnitude of the average local superfluid density if we assume that 
only a depends on T. We remark that the possibility of a contribution to the entropy of 
this type is discarded from the beginning in a frustrated XY or London model description 
of the phase transition@iH, since the magnitude of the order parameter is implicitly held 
constant and no attempt is normally made to estimate the temperature dependence of the 
parameters of this effective model. 



Using Eqs. (Wm), (H) and fl24|) , the latent heat can be expressed in terms of material 



parameters and discontinuities in intensive thermodynamic variables across the transition. 

T m 5S = N,N z k B T m g 2 m [- ° < 0> +%6<Ej> + nMF ™ „ " 6 < A >]. (25) 

Z Z JrL c2 {l m ) — n m al 

(Recall that Ej = AoT.) The magnetization in the LLL model comes entirely from the 
dependence of the model parameters on field and is proportional to Ao: 

AnM _kj£l_ Ap _ (H% F - H)A 

H 2ir£ 2 d H*£ F - H Ak 2 H ' 1 ' 

Since A increases on freezing, the internal magnetic field decreases and the vortex lattice 
therefore expands slightly. This expansion will always occur in the LLL limit @ Eq. fl2"6| ) 
is valid for extreme type-II superconductors in the LLL regime since the screening current 
gives rise to a small correction to the internal field which is approximately uniform spatially. 
There is no need to distinguish between internal and external fields on the right hand side of 
this equation. (This is the same property which has allowed us to ignore thermal fluctuations 
in the internal magnetic field.) 

Our calculations in principle allow theory and experiment to be quantitatively compared 
for any strongly anisotropic layered superconductor. For a specific superconducting mate- 
rial, with known H^ F (T), k, 7, d, and do, the point {g m ,rj m ) in parameter space where 
an apparent first order phase transition occurs can be calculated from field strength and 
temperature using Eq. (|TJ) and Eq. (|T3[) . The applicability of the LLL model can be judged 
from the proximity of the melting point to the estimates given here in Fig. Our estimate 
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for the latent heat and the magnetization discontinuity can then be calculated using Eq. fl25|) 
and Eq. (|26| ) and the discontinuities in Ej and A plotted in Fig. [7]. We remark that where 
the LLL limit is valid, the relative change in the internal magnetic field for strong type-II 
superconductors at the transition is always small, less than ~ 1/(200k 2 ). The entropy jump 
per vortex per layer can be substantially larger and may be estimated to be: 

ss gls {Ap) 

kBN+N, ~ (1 - T c0 /T m )(l - H/HM F (T m )) 1 ] 

The entropy jump per vortex, per layer can be much larger than the value estimated in XY 
model simulations E3 in which the dominant contribution to Eq. (|25|) is absent. We compare 
our results with some recent experiments in a subsequent subsection. 



C. Off-diagonal long-range-order 

We now turn our attention to the evidence for superconductivity and off-diagonal or- 
der that comes from simulations. In Fig. |8] we compare results for local superfluid density 
correlation functions and off-diagonal correlation functions at the phase transition point 
where contributions due to both ordered and disordered states contribute. The peak in the 

n+4 (g = 0) distribution function^ centered on a non-zero value comes from the same 
ordered state configurations that give peaks in the probability distribution functions for 
S n ,n(G, 0) and S ntn+4: (G, 0). Note that phase fluctuations along the field direction cause this 
quantity to have a non-zero imaginary part, and that the peak in the distribution function 
has a greater width in the imaginary direction than in the real direction, indicating that 
phase fluctuations have a greater importance than amplitude fluctuations. There is no ev- 
idence of these off diagonal correlations above the vortex lattice melting transition. This 
result contrasts with XY model studies which find a superconducting vortex line liquid 
state that has phase coherence along the field direction without broken translational sym- 
metrySEI. The explanation for this difference might originate in the fact that no clearly 
defined region which is outside of the vortex cores exists in the LLL approximation, mak- 
ing liquid state vortex position fluctuations more effective in destroying coherence along the 
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field. The absence of such a phase in the LLL approximation does not rule out the possibility 
that it exists at weaker magnetic fields. 

The concept of ODLRO and has long played a fundamental role in the theory of super- 
fluidity and superconductivity.^ Only recently^E however, have attempts been made to 
apply this concept to superconductors in a magnetic field. These efforts have led to some 
confusion and controversy, partly due to the gauge dependence of Cooper pair correlation 
functions in a magnetic field. The relationship is an interesting one in GL models generally 
and in LLL models in particular.il (The relationship is partially obscured in XY and Lon- 
don models for vortex states.) In the absence of a field, superconductivity in 3D systems 
is associated with Bose condensation of Cooper pairs, i.e. with macroscopic occupation of 
their zero momentum state. This properties immediately implies ODLRO , i.e. that 

lim (*(r + R)V(r)) ^ 0. (28) 

R^oo 

In a magnetic field, Cooper pair states are labeled by a momentum along the field direction 
(restricted to the appropriate Brillouin zone for the Lawrence-Doniach model), by a Landau 
level index, and by a label that counts the states within a Landau level. Superconduc- 
tivity is presumably associated with the macroscopic occupation of a Cooper pair state with 
k z = and with the lowest Landau level index, but it less obvious what state should be 
occupied within a Landau level. 

In order to avoid the requirement for prior knowledge about the state in which Bose 
condensation has occurred, we propose defining off-diagonal long-range-order in a magnetic 
field in terms of the eigenvalues, p a , of the Cooper pair density matrix. For the case of the 
Lawrence-Doniach model the density matrix takes the form 

(*„(f)*»(*0> = Y.PMr'W n {^) (29) 

a 

where (fi a is a normalized eigenfunction of the density matrix. (J2 n J d 2 \4>n(^)\ 2 = !•) ^aPa 
will be proportional to the system size whenever superconducting fluctuations are present 
in the system. For example, in our finite-size LLL Lawrence-Doniach simulations, 
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2^Pa = -p A . (30) 

In a magnetic field, a GL model may be said to possess ODLRO when at least one eigen- 
value of the density matrix is extensive and the corresponding eigenfunction is extended. 
We remark that in a magnetic field this criterion cannot be satisfied unless the system 
breaks translational symmetryU either spontaneously, or as a consequence of pinning inho- 
mogeneities. Convincing arguments can be advanced^ in support of the view that super- 
conductivity, both parallel and perpendicular to the field, will occur i f§ the system has this 
ODLRO property. 

In Abrikosov's mean-field theory, translational symmetry is broken to form a triangular 
vortex lattice, and the density matrix has a single extensive non-zero eigenvalue. When 
fluctuations are included ODLRO will survive if one eigenvalue of the density matrix remains 
macroscopically larger than all others. This quantity is the analog for GL models of the 
condensate fraction in a boson particle superfluid. In Fig. ||]we plot Z, defined as the ratio of 
the largest eigenvalue of the density matrix to the sum of all eigenvalues, as a function of g 
at rig = —0.045. For this value of rig, a first order melting transition occurs at g ~ — V20.52. 
The insets in this figure show the systems size dependence of the 'condensate fraction' for 
temperatures on the vortex solid and vortex liquid sides of the phase transition. From 
Fig. |^ we see that our calculations indicate that ODLRO occurs simultaneously with broken 
translational symmetry, at least when the LLL approximation is valid. 

D. Comparing with experiment 

Since the discovery of high temperature superconductivity, the superconducting transi- 
tion in a magnetic field has received considerable experimental attention. It is abundantly 
that the mean field phase transition, particularly as seen in specific heat and magne- 
tization measurements, is spread over a wide fluctuation regime and looks broadly similar 
to what would be expected if the system were below its lower critical dimension and no true 
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phase transition occurred. However, according to XY model simulationsU'El and this and 
other simulations of the LLL model, there is a weak first order phase transition which occurs 
well below the mean-field transition temperature and is accompanied by a relatively small 
latent heat and magnetization discontinuity. In this section we compare some of the more 
recent experimental studies!^ of systems that are free of pinning inhomogeneities with our 
simulations. In making this comparison we are seeking quantitative agreement which would 
indicate that the LLL model provides a good description of the experiment in question. 

We first compare with experimental studiesiH'iS of YBCO which have identified a first 
order phase transition in the field and temperature range (H m ,T m ) = (0-6 Tesla, 93-81 K) 
revealeoffl by a small reversible discontinuous increase in the diamagnetism of the system 
with decreasing temperature. The latent heat of the transition can then be determined from 
the temperature dependence of the critical field using the appropriate Clapeyron equation. 
Between 2 Tesla and 5 Tesla, the magnetization jump increases from ~ 0.012 Gauss to 
~ 0.025 Gauss so that the relative decrease in the internal magnetic field, 4nAM/H, varies 
slowly and is close to 7 x 10~ 6 . Over the same field range, the latent heat per vortex per 
double-layer of YBCO is fairly constant and close to 0.65fcsT m . From Table 1 we see that 
in YBCO gr] is ~ —0.3 at H = 5 Tesla and varies as if -1 in the temperature interval of 
interest. Thus the highest fields studied in the experiment are slightly outside the portion of 
phase diagram we explored here. Nevertheless, the experimental melting phase diagram for 
fields larger than 2 Tesla is consistent with the extrapolation of our phase diagram to larger 
r\. This claim is reinforced by the LLL Lawrence-Doniach simulations of Sasik and Stroud 
which specifically modeled YBCO and found excellent agreement with the measured phase 
boundaryU Using material parameters from Table 1 for YBCO, dH^ F /dT = —2.0 Tesla/K 
and the jump in A in Fig. [7| gives 4nAM/H ~ 8.9 x 10~ 6 at H = 5 Tesla. Similarly the 
latent heat calculated from Eq. (^) is 0.58k B T m per-vortex per double-layer; both values 
agree well with experiment. For YBCO our results confirm the excellent agreement with 
experiment found by Sasik and Stroudil for the magnetization discontinuity. However, these 
authors seriously underestimated the latent heat by neglecting the last term in Eq. ( pq) which 
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contributes ~ 90% of the total value. We remark that the latent heat comes primarily from 
the change in entropy content at microscopic length scales associated with a change in the 
magnitude of the superconducting order parameter and not from changes in the entropy 
content of vortex configurations. In the light of this remarkable quantitative agreement 
between theory and experiment, it seems clear that above H = 2 Tesla the phase transition 
cannot be reliably described by London or XY models. 

We now turn to the experimental studies! of the much more anisotropic BSCCO system 
where first order phase transitions have been found in a much weaker field range: (H m , T m ) = 
(0-380 Gauss, 90-40 K). The latent heat at the transition is again ~ LOfceT™ per vortex 
per layer and the jump in the internal magnetic field (4ttAM) on freezing is ~ 0.3G. Con- 
verting (H m ,T m ) to our dimensionless coupling parameters g m and r] m , we find g m values 
up to ~ 100 for BSCCO in the ranges of fields probed experimentally. This is clearly in- 
consistent with our LLL numerical results for the melting temperature. We have concluded 
that in BSCCO, weak interlayer coupling contributes to the high transition temperatures 
and relatively low superfluid densities present in all high T c superconductors, resulting in 
fluctuations that are so strong that the melting transition occurs outside the validity range 
of the LLL model, even at temperatures quite far from the zero-field transition temperature. 
In accepting this conclusion, it is important to realize that the largest part of the entropy 
reduction associated with ordering appears in the specific heat and in the magnetization 
at temperatures centered on the mean-field transition temperature. At moderate fields, it 
has already been convincinglyS demonstrated that thermodynamic properties across these 
broad fluctuation regimes are well described by the LLL model; the melting transition occurs 
near the low temperature extreme of the fluctuation regime and it is only here that the LLL 
approximation fails for BSCCO and similar strongly anisotropic systems. 
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V. SUMMARY 



We have presented a discussion of the physics of vortex states of layered superconductors 
based on Monte Carlo simulations of the lowest Landau level limit of the Lawrence-Doniach 
model. The Lawrence-Doniach model provides a realistic description of these systems and 
may be used to confront theory and experiment on a quantitative level. The lowest Landau 
level limit of the model is valid at temperatures close to the mean field transition and at 
stronger fields. In this limit we find that the model has a single first order phase transition 
between a high temperature phase with no broken symmetries and a low-temperature phase 
with broken translational symmetry and off-diagonal long-range-order, which we discuss in 
terms of the eigenvalue distribution of the Cooper pair density matrix. 

Although it has been difficult to accurately delimit the validity of the lowest Landau 
level constraint on purely theoretical grounds, quantitative comparison with experiment 
can provide a pragmatic criterion. With this motivation, we have determined a universal 
phase diagram for this limit of the Lawrence-Doniach model in terms of two dimensionless 
parameters which can be converted into field and temperature using the material param- 
eters of a specific layered superconductor. Our numerical simulations have focused on the 
strongly anisotropic limit where the discreteness of the superconducting layers plays a role. 
Comparison with experiments indicates that the LLL approximation is accurate along the 
phase boundary for fields above ~ 2 Tesla in YBCO, but that it is not accurate at more 
strongly suppressed transition fields in BSCCO systems. In YBCO we find good agreement 
between experiment and simulations for both magnetization jumps and latent heats associ- 
ated with the first order phase transition. We conclude that the lowest Landau level limit 
of the Lawrence-Doniach model provides a quantitatively reliable description of the phase 
transitions studied experimentally in these systems. In this model the latent heat at the 
phase transition is primarily due to entropy changes on microscopic length scales which are 
captured by the temperature dependence of the model parameters and are not accurately 
modeled in XY or London descriptions. 
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TABLES 

TABLE I. Material parameters and corresponding gr] at typical temperatures and magnetic 
fields for several superconductor systems with strong planar anisotropy. 



Material 


K 


7 


d (A) 


d(A) 


K 1 2 d 2 d Q 112 


gv 


BSCCoi 


100 


55 


4 


15 


3.40 x 10 7 


0.0075 (at H = 1 T, T = 50.0 K) 


YBCO0 


65 


7.7 


3 


11.4 


2.89 x 10 5 


0.30 (at H = 5 T, T = 82.5 K) 


MoGell 


140 


22.4 


60 


125 


1.41 x 10 8 


0.013 (at H = 1 T, T = 1 K) 


NbGe@ 


83.4 


4.6 


180 


202 


5.37 x 10 6 


0.33 (at H = 1 T, T = 1 K) 



30 



FIGURES 



0.05 




0.00 2.00 4.00 6.00 8.00 10.00 



9j 

0.80 i - 




0.0 2.0 4.0 6.0 8.0 10.0 



9j 

FIG. 1. (S n ' %n (q x , 0)) at i] = 0.01 and g = —t/10 (upper panel), and g = —a/30- These 
results are for a finite system with N$ = 36 and N z = 8 and were obtained by averaging over 10 6 
Monte-Carlo steps. 
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FIG. 2. Landau-Ginzburg energy distribution function at the finite system phase transition 
point for various system sizes. Energies are in units of the mean-field condensation energy, 
N^N z k B Tg 2 /(3 A . 
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FIG. 3. Two dimensional contour plot for the distribution function of (E2d,Ej) at the finite 
system phase transition point (rj m ,g m ) = (20.52,-0.01). The peak located at smaller < E^d > 
and < Ej > comes from 3D vortex solid configurations. The peak located at larger < E^d > and 
< Ej > comes from vortex liquid configurations. 
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FIG. 4. The phase diagram rj m = f(g m ) in terms of the two dimensionless parameters g and 
T) we use to characterize intra-layer and inter-layer couplings. Phase transition points located in 
our simulations are indicated by solid circles. This phase boundary was obtained by spline fitting 
using the four smallest 77 phase transition points. The two points at larger 77 are ones at which 
the system appears to undergo a phase transition as judged by the temperature dependence of 
the correlation functions, but we are not able to locate the double peak structures due to limited 
computation time and small finite system size. 
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FIG. 5. Landau-Ginzburg energy distribution function near the finite system phase transition 
point for g = —9.5 and rj = 0.06. These results are for N z = 32 and N<p = 16. Energies are in units 
of the mean-field condensation energy, N<pN z kBTg 2 / '(3a- 
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FIG. 6. The g dependence of thermodynamic variables Aq, Pa and T for gr\ = —0.01. The 
calculations are done for the system size iV^ = 16, N z = 8. 



36 




37 



1.5 



_ 1-0 

o 

(3 

CO 

^ 0.5 




0.0 0.2 0.4 0.6 0.8 1.0 

S„,„(G,0) 





0.2 






0.1 












p 








+ 


0.0 


\\ 


n ■ 


| 


-0.1 








-0.2 





-0.5 0.0 0.5 1.0 

Re(S„, n+4 (C,0) 

1 .0 i - 




-1 o r . .................. 1 

-1.0 -0.5 1 0.0 0.5 1.0 
R e(A n , n+4 (0)/A ) 

FIG. 8. The top panel shows the distribution function for S n:Tl (G,0); the middle panel is a 
contour plot of the distribution function for 5„ jr! ,+4(G, 0) in the complex plane; the bottom panel 
is a contour plot of the distribution function for A njn+ 4(cf = 0)/Ao in the complex plane. All plots 
are at the melting transition point of a finite system with iV = 6x6x8, g = — V20.5, and r/ = 0.01. 
G is a member of the first shell of reciprocal lattice vectors of the Abrikosov lattice. 
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FIG. 9. The ratio of the largest eigenvalue of the density matrix to the sum of all eigenvalues, 
as a function of g at rjg = —0.045. For this value of grj a first order melting transition appears 



to occur at g ~ —y/20.5. The size dependence of Z at particular temperatures on the vortex solid 
and vortex liquid sides of the phase transition are shown in the right and left insets respectively. 
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